The thermodynamic hypothesis proposed by Anfinsen is that the native state tends to have the lowest free energy.^[42] To evaluate an RNA 3D structure, the proposed potential considers two energy terms. The total energy u_total of a structure is given by

The energies between the two atoms in Eq. (1) are knowledge-based potentials derived from a Boltzmann distribution or Bayesian formulations,^[43] and the energy between two atoms with type i and type j is

The parameters of the proposed potential were calculated based on the statistical analysis of known RNA 3D structures in the PDB/NDB database. First, we gathered 1369 non-redundant structures including most of the common RNAs (e.g., mRNA, rRNA, tRNA, ribozyme, and riboswitch) and various RNA complexes where only RNAs were taken into account. Second, we discarded the structures with sequence identities greater than 80% using the BLASTN program with default options^[44] and removed low-quality structures with a resolution . Finally, 380 structures with only normal heavy atoms containing various typical motifs (e.g., hairpin/internal/bulge loops, double/triple helices, junctions, and non-Waston–Crick base pairs) were selected as the training set and the PDB ID for each of selected structures is listed in Table S1. Based on the careful selection of the training set, the distance-dependent parameters of the potential can be obtained as a matrix with a bin width of 0.3 Å and a cut-off distance of 20 Å as used in similar recent potentials.^[30,33]

In structure ranking, the initial potential obtained above can be further optimized for decoys of an RNA by a retraining mechanism (Fig. S5). In the retraining process of each RNA, based on the energies calculated by the initial potential, structures with the lowest energies are used as a new training set to further obtain an RNA-specific potential, and energies of all conformations can be calculated again by the retrained potential (Fig. S5(a)). Note that the number of conformations with low energies is generally taken as 10, where the potential has the highest precision (Fig. S5(b)).

To test the performance of the proposed potential, we used two different decoy sets widely used in recent works to analyze the performance of knowledge-based potential.^[8,29–33] Test set I, which has been used by RASP and 3dRNAscore, consisted of decoys generated from 85 native structures (only with normal heavy atoms) by MODELLER^[45] with a set of Gaussian restraints for dihedral angles and bond stretches and can be downloaded from http://melolab.org/supmat/RNApot/Home.html. Test set II contained three parts which were a set of decoys generated by three separate methods recently used to assess the scoring functions of KB, Rosetta, and 3dRNAscore. The first part of Test set II contained five distinct RNA decoy sets generated by MD by using a position restraint potential on each heavy atom to constrain the motions of the RNA, named MD decoys.^[46,47] The second part of Test set II included decoys of 15 RNAs generated by the normal-mode (NM) perturbation method,^[48] named NM decoys. These two decoys can be downloaded from http://csb.stanford.edu/rna/. The third part of Test set II, named FARNA decoys, included 3D structures of 17 RNAs with many non-native base pairs blindly predicted by Rosetta 3.1 from sequences downloadable from https://daslab.stanfold.edu/resources. It should be noted that a very small portion of tested RNAs existed in the training set (Sup. Table 1) that could not affect the reliability of the proposed potential.

Generally, the enrichment score (ES) is often used to describe the performance of a scoring function in identifying the best RNA structures.^[32,49] The ES is defined as

Owing to the different sample sizes for the two potentials in Eq. (1), the coupling parameter ε is extremely important in the combination of the two energy terms.^[50] To obtain a proper value of ε, we selected five RNAs (PDB ID: 1KKA, 1MSY, 1ZIH, 255D, and 437D) with different sizes and representative motifs including hairpins, bulges, pseudoknots, and junctions. For each of the five RNAs, we used the MC-Fold/MC-Sym pipeline to generate a series of 3D structures^[7] that were scored by the proposed potential with different values of ε. As shown in Fig. S2, when ε is equal to 0.15, the calculated average ES (RMSD) of the five decoy sets was the maximum value. Therefore, an ε value of 0.15 was used for all the tested RNAs in this study.

To test the performance of the proposed potential, we first use our potential to select the native structure from decoys based on the thermodynamic hypothesis.^[42] For decoys in Test sets I and II and the corresponding native structures, we use the proposed potential with initial parameters to calculate the energy for each conformation and select the conformation with the lowest energy as a native structure. As shown in Fig. 1 and Figs. S3 and S4, the proposed potential accurately selected the native structures from even the poor decoys. For example, our potential identified 83 native structures out of 85 RNAs in Test set I and 34 native structures from 39 RNAs in Test set II.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) Counts of native states identified correctly from Test sets I and II by RASP, KB, Rosetta, and our proposed potential.

We also compared the proposed potential against well-established existing potentials such as RASP, KB, and Rosetta in identifying the native structures from different decoys. We employed the RASP, KB, and Rosetta potentials with corresponding programs or software to select native structures for the two test sets, and the results are shown in Fig. 1. For Test set I, 77, 80, and 53 out of 85 native structures were identified by the three potentials, respectively, which are all less than that identified by our potential. Although the proposed potential was as good as the RASP for Test set II in selecting native structures, it performed slightly better than the other two potentials.

It is noted that several native structures, such as 1ESY, 1KKA, and 1QWA in the test sets, could not be correctly selected by either the proposed potential or the other existing potentials as shown in Fig. 1 and Figs. S3 and S4. One possible reason could be that these structures have their own special features, such as complex ligand binding or structure stretching.^[51–53] In addition, the number of RNAs in the training set was limited and resulted in over-training.^[34]

The identification of the best structures is a vital task in RNA tertiary structure prediction.^[8,13] For example, although a recent new method using evolutionary restraints of RNAs was proposed to predict RNA 3D structures, it still needs to recognize the best conformations from the predicted structures using a knowledge-based potential such as 3dRNAscore.^[21] However, it is still challenging and essential for the existing predictive methods to accurately recognize the best conformations from a series of predicted conformations. Thus, the precise identification of the best RNA conformations from a set of near-native structures is an essential purpose for the scoring function. To assess the quality of our potential, we employed the proposed potential to discriminate and rank the near-native conformations of RNA.

For each RNA in Test set II, we first calculated the energies of all decoys excluding the native structure by using the proposed potential^[30,33] and ranked the decoys based on energy (Fig. S5). Then, selecting 10 decoys with the lowest energies as a new training set, the potential was trained again to consider the characteristics of the RNA^[37,39] This retraining process step proved to be effective for improving the accuracy of RNA structure evaluation. Finally, we further ranked the decoys of the RNA again based on the energies scored by the retrained potential (Figs. S6 and S7) and assumed that the lower the energy, the better the structure.

As shown in Fig. 2 and Figs. S6 and S7, for each set of decoys from different RNAs, the energy scores generally declined with a decrease of RMSD (or DI) and such a reduction became steeper near small values of RMSD, which suggests that the proposed potential was very effective at ranking near-native structures and selecting the best.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) Energy as a function of RMSD for five typical MD/NM decoys (rows). The proposed potential, 3dRNAscore, KB, RASP, and Rosetta are shown in turns (columns).

In addition, we benchmarked the performance of 3dRNAscore, RASP, KB, Rosetta, and the proposed potential in selecting the best structures on Test set II using both RMSD-based ES and DI-based ES. As shown in Fig. 3 and Table S2, the overall results of Test set II by the proposed potential (ES (RMSD) = 5.6 and ES (DI) = 5.5) were better than those from the KB (ES (RMSD) = 3.7 and ES (DI) = 3.7), RASP (ES (RMSD) = 3.6, and ES (DI) = 3.8), 3dRNAscore (ES (RMSD) = 4.5 and ES (DI) = 4.6), and Rosetta (ES (RMSD) = 2.7 and ES (DI) = 2.7), which is also suggested by the p-values less than 0.01 from the non-parametric and distribution-free Kolmogorov–Smirnov (KS) test (Table S3(a) and S3(b)).

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) Comparison between average values of RMSD-based ES (a) or DI-based ES (b) calculated from Rosetta, 3dRNAscore, RASP, KB, and our potential with retraining mechanism on Test set II: all decoys, FARNA decoys, NM decoys, and MD decoys.

Figure 3 and Table S2 also show that for FARNA decoys, the ES values calculated by the proposed potential (ES (RMSD) = 2.7 and ES (DI) = 2.6) as well as the others (ES ) were smaller than those of the other decoys. This may be attributed to the fact that the FARNA decoy set contained many tertiary structures that included many non-canonical base pairs.^[8] The result also suggests that all the existing potentials’ evaluation of structures predicted from blind prediction experiments (e.g., RNA-Puzzles) needs to be improved.^[5] In future studies, the proposed potential should be further tested in more real RNA prediction cases.

Generally, base-pairing and base-stacking interactions provide a strong force in stabilizing RNA 3D structures.^[17,54] Thus, the capacity to precisely capture base-pairing/stacking interactions is one of the significant criteria for determining the quality of potentials. We further analyzed the energies between the specific atom pairs involved in base-pairing and base-stacking such as N1 of adenine and N3 of uracil in the proposed potential (Fig. 4(a)) which were deduced from a Boltzmann distribution or Bayesian formulations (Eq. (1)) based on our training set.^[42,43] As shown in Fig. 4(a), there were two visible wells: one at 3.0 Å and the other at 5.1 Å, which accurately represent the distances between corresponding atoms in pairing and stacking bases, respectively,^[55,56] as shown in Fig. 4(b). We also calculated the energies between any two bases using the proposed potential for RNA (PDB ID: 434D) because RNA has a typical helical structure and represents many RNAs.^[8,57] The secondary and 3D structure of RNA are shown in Fig. 4(d). As shown in Fig. 4(c), the lower the energy, the stronger the base-base interaction: (i) energies between the nucleotides 1–14, 2–13, 3–12, 4–11, 6–9, and 5–6 are the lowest energies, which exactly represents the four types of Watson–Crick base pairs in 434D (Fig. 4(d)), and (ii) energies between bases 3–12, 4–11, 6–9, and 5–6 are stronger than those between bases 1–14 and 2–13, which represents an interaction of cytosine and guanine stronger than that of adenine and uracil and is in line with previous works.^[8,57] Above all, the proposed potential effectively captured the interaction of base-pairing and base-stacking.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) (a) Energy distribution of distance between N1 of adenine and N3 of uracil. (b) Diagram of representative distance between N1 of adenine and N3 of uracil. (c) Base-stacking and base-pairing energies between each possible base-pair calculated by the proposed potential for 434D. (d) Secondary structure and 3D structure of 434D.

Based on the carefully selected training set shown in Figs. 1–3, the proposed potential with the combination of the new energy term u_intra and the employment of a retraining mechanism performed better than the existing methods in assessing RNA structures for Test sets I and II.

To clarify the contributions of the two improvements, we further performed two additional tests for each RNA decoy in Test set II using our potential without the u_intra or the retraining mechanism (Table sup). As shown in Fig. 5, the two improvements made positive contributions to the overall improvements, which is also suggested by the p-values less than 0.01 from the non-parametric and distribution-free KS test (Table S5), and the inclusion of the retraining mechanism made a very strong contribution (Fig. 5). As shown in Table S5, the proposed potential with only U_inter performed slightly better than the other potentials, which possibly contributed to the success of our training set in including extensive non-redundant RNAs and embracing most motifs. It is noted that although the local interactions U_intra introduced into the proposed potential only slightly improved the results (Table S5), it would be more significant in the evaluation of complex RNAs including many local tertiary interactions.^[33] Although the addition of the new term slightly improved the performance of our potential in some cases (Table S5), it is still particularly necessary for large RNAs with complex structures because the new term represents the local geometrical features and flexibility of RNA structures^[40,41] Overall, the specificity of RNA molecules and the interactions between two atoms within one nucleotide can be necessary.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) Comparison between average values of RMSD-based ES calculated from the proposed potential, the potential without uintra, and the potential without retraining mechanism on Test set II.

However, for some RNAs such as 1DQF, 1I9X, 1J6S, 1KD5, 1KKA, and 1MHK in the FARNA decoy set, the employment of the retraining process can result in a decrease of ES (Table S4). The possible reason is that the retraining process is based on the results generated from the initial scoring; however, for these RNAs with complex tertiary interactions (e.g., non-native hydrogen bonds), the initial scoring obtained by the proposed potential as well as the other existing potentials still needs improvement.

Although the knowledge-based statistical potential has proven to be an efficient method for structural evaluation, the selection of an effective reference state is still an inherent limitation for the knowledge-based scoring function. Selecting diverse reference states to score the same decoys can lead to different results.^[58,59] Since an ideal reference state is not achievable, current statistical potentials generally construct reference states through randomizing disconnected atoms; e.g., the averaged reference state in 3dRNAscore and RASP and the quasi-chemical approximation reference state in KB, which ignore the diversity of various RNAs.^[60] Therefore, on the basis of choosing the averaged reference state, we further propose a retraining mechanism in the proposed potential to distinguish the characteristic of different RNAs, which has proved to be very effective in assessing RNA 3D structures (Fig. 5).